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ABSTRACT 

Ill-posed linear inverse problems (ILIP), such as restoration and re- 
construction, are a core topic of signal/image processing. A standard 
approach to deal with ILIP uses a constrained optimization prob- 
lem, where a regularization function is minimized under the con- 
straint that the solution explains the observations sufficiently well. 
The regularizer and constraint are usually convex; however, sev- 
eral particular features of these problems (huge dimensionality, non- 
smoothness) preclude the use of off-the-shelf optimization tools and 
have stimulated much research. In this paper, we propose a new effi- 
cient algorithm to handle one class of constrained problems (known 
as basis pursuit denoising) tailored to image recovery applications. 
The proposed algorithm, which belongs to the category of augmented 
Lagrangian methods, can be used to deal with a variety of imaging 
ILIP, including deconvolution and reconstruction from compressive 
observations (such as MRI). Experiments testify for the effectiveness 
of the proposed method. 

Index Terms — Optimization, inverse problems, image recon- 
struction/restoration, compressive sensing, total variation, tight frames. 

1. INTRODUCTION 
1.1. Problem Formulation 

Linear inverse problems constitute one of the central themes of sig- 
nal/image processing. In this class of problems, a noisy indirect ob- 
servation y, of an original signal x, is modeled as 

y = Bx + n, (1) 

where B is the matrix representation of the direct operator and n 
is noise. In the sequel, we denote by n the dimension of x, thus 
x £ R n , while y G R m . In the classical problem of image deblur- 
ring/deconvolution, B is the matrix representation of a convolution 
operator. In other reconstruction problems, B represents some linear 
direct operator, such as of tomographic projections (Radon trans- 
form) or a partially observed (e.g., Fourier) transform (as in com- 
pressive MRI fl9l ). 

Usually, the problem of estimating x from y is ill-posed (e.g., if 
m < n), thus requiring some sort of regularization. In the presence 
of noise, a natural criterion to infer x from y has the form 171 1201 

min<j!>(x) subjectto ||Bx — y||a < £, (2) 

where <p : R" — *• R is the the regularizer and e > a parameter 
which depends on the noise variance. In the case where <f>(x) = 
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||x|| i, the above problem is usually known as basis pursuit denoising 
(BPD) |8|. The basis pursuit (BP) problem is the particular case of 
l|2j for e = 0. In recent years, an explosion of interest in problems of 
the form ([2} was sparked by the emergence of compressive sensing 
(CS) [5], [ 10 1 . The theory of CS provides conditions (on matrix B 
and the degree of sparseness of the original x) under which a solution 
of (2) is an optimal (in some sense) approximation to the "true" x. 

In most image recovery and CS problems, the regularizer <j> is 
convex but non-smooth; typical examples are the total variation (TV) 
[5|, |23 | and^i norms. Problem |2]( is thus convex, but the very high 
dimension (usually > 10 4 ) of x and y precludes the direct applica- 
tion of off-the-shelf optimization algorithms. This difficulty is fur- 
ther amplified by the fact that matrix B only "exists" as an operator; 
i.e., there are efficient algorithms to compute products of B (or B T ) 
by some vector (image), but it is highly impractical to extract and 
manipulate individual blocks, rows, or columns of this matrix. 

1.2. Previous Work 

Most state-of-the-art methods for dealing with linear inverse prob- 
lems, under convex, non-smooth regularizers (namely, TV and 
consider, rather than |2j, the unconstrained problem 

min -||Bx - y\\\ + r 0(x), (3) 
x 2 

where r 6 R+ is the so-called regularization parameter. Of course, 
problems {2]l and ([3j are equivalent, in the following sense: for any 
e such that problem |2| is feasible, a solution of Q is either the null 
vector, or else is a solution of ([3}, for some r 1 15 1. 

The currently fastest (publicly available) algorithms for solving 
p), include: gradient projection for sparse reconstruction (GPSR) 
1151 ; fast iterative shrinkage/thresolding algorithm (FISTA) [ Q; two- 
step 1ST (TwIST) (2); and sparse reconstruction by separable ap- 
proximation (SpaRSA) |28|. These methods were shown to be con- 
siderably faster than earlier methods, including ll_ls 1 18] and the 
codes in th e l\-magic (http : / / www . 11 -magic ■ org) and the 
SparseLab (http : / /sparselab . Stanford, edul toolboxes. 
Very recently, we have introduced a new algorithm, called SALSA 
(split augmented Lagrangian shrinkage algorithm); experiments on 
a set of standard image recovery problems show that SALSA is faster 
than GPSR, TwIST, FISTA, and SpaRSA fl3l . 

Although it is usually easier/simpler to solve an unconstrained 
problem than a constrained one, formulation |2| has an important 
advantage: parameter e has a clear meaning (it is proportional to the 
noise variance) and is much easier to set than parameter r in {3}. Of 
course, one may solve |2]( by using one of the algorithms mentioned 
in the previous paragraph to solve j3} and searching for the "correct" 



value of r that makes {3} equivalent to |2j. Clearly, this is not effi- 
cient, as it involves solving many instances of d3). Obtaining fast 
algorithms for solving |2| is thus an important research front. 

There are few efficient algorithms to solve {2} in an image re- 
covery context: x and y of dimension > 10 4 (often > 10 6 ), B 
representing an operator, and (f> a convex, non-smooth function. A 
notable exception is the recent SPGL1 |26|, which (as its name im- 
plies) is specifically designed for i\ regularization (0(x) = ||x||i). 
Other methods for solving problems with the form l|2j, for equal 
to the l\ or TV norms, are available in the t\ -magic package; how- 
ever, as shown in |26|, those methods are quite inefficient for large 
problems. General purpose methods, such as those in the SeDuMi 
package jhttp : / / sedumi . ie . lehigh . edu| l, are simply not 
applicable when B is not an actual matrix, but an operator. 

The Bregman iterative algorithm (BIA) was recently proposed 
to solve |2| with e = 0, but is not directly applicable when e > 
1291 . To deal with the case of e > 0, it was suggested that the 
BIA for e = is used and stopped when ||Bx — y || 2 < e (4|, |29|. 
Clearly, that approach is not guaranteed to find a good solution, since 
it depends strongly on the initialization; e.g., if the algorithm starts 
at a feasible point, it will immediately stop, although the point may 
be far from a minimizer of <j>. 

1.3. Proposed Approach 

In this paper, we introduce an algorithm for solving optimization 
problems of the form |2j. The basic ingredients are the following: 
the original constrained problem (|2| is transformed into an uncon- 
strained one by using an indicator function of the feasible set; the 
resulting unconstrained problem is transformed into a different con- 
strained problem, by the application of a variable splitting opera- 
tion; finally, the obtained constrained problem is attacked with an 
augmented Lagrangian (AL) scheme [22], which is a variant of the 
alternating direction method of multipliers (ADMM) II II . Since (as 
SALSA), the proposed method uses variable splitting and AL opti- 
mization, we call it C-SALSA (for constrained-SALSA). 

The resulting algorithm is more general than SPGL1, in the 
sense that it can be used with any convex regularizer (f> for which 
the corresponding Moreau proximity operator |9|, defined as 



1 2 
*r^(y) = argmin-||x-y|| 2 + T(i 
x 2 



(4) 



has closed form or can be efficiently computed. Below, we will 
show examples of C-SALSA where x is an image, <f> is the TV norm 
1231 . and is computed using Chambolle's algorithm (6). An- 
other classical choice is </>(x) — ||x||i, which leads to ^ r T ^>(y) = 
soft(y, r), where soft(-, r) denotes the component-wise application 
of the soft-threshold function y t—> sign(y) max{jj/j — r, 0}. 

C-SALSA is experimentally shown to efficiently solve image re- 
covery problems of the form {2|, such as MRI reconstruction from 
CS-type partial Fourier observations using TV regularization. More- 
over, C-SALSA is also shown to be faster than SPGL1 in wavelet- 
based image deconvolution problems under l\ regularization. 

The paper is organized as follows. Section[2]brieny reviews vari- 
able splitting and ADMM. Section[3]contains the derivation leading 
to C-SALSA. Section [4] reports experimental results, and Section [5] 
ends the paper with a few remarks and pointers to future work. 

2. VARIABLE SPLITTING AND ADMM 



Consider an unconstrained optimization problem 
min/ 1 (u) + / 2 (Gu), 



where G G R xn . Variable splitting (VS) is a simple procedure 
that consists in creating new variables, say v and w, to serve as the 
argument of each of the terms, fi and / 2 , under the constraints that 
w = u and v = Gu, that is, 



min /i (w) 

u,wgK",v6M d 



Mv) 



subject to 



w = u 
v = Gu. 



(6) 



Problem ([6j is clearly equivalent to the unconstrained problem {3J. 
The rationale behind VS is that it may be easier to solve the con- 
strained problem |6]) than to solve its unconstrained counterpart |5j. 
It is important to stress that the VS in l[6} is not the one commonly 
used, where only variable v is created; however, as shown below, the 
proposed VS will lead to a very effective algorithm. 

Other variants of VS were recently used in several image pro- 
cessing problems: in |27|, it was used to obtain a fast TV-based 
algorithm; in |3 |, it was used to handle problems with compound 
regularization. VS also underlies the recent split-Bregman methods 
1161 . but there the splitting is different and with a different goal. 

Using an augmented Lagrangian (AL) approach to handle prob- 
lem |6]( leads to the following algorithm, also known as the method 
of multipliers (MM) Q7), |24) (see also 02), for details): 

(Uk+i, Vfc+i, Wjt+i) G arg min \ /i(w) + /a(v) + 

|Gu-v-b fc |||' + ^||u-w-c fc |||} (7) 



2 



b fe+ i = b fc + Gu fe+ i - Vfc+i (8) 
Cfc+i = c*; + Ufe+i - Wh+i- (9) 

Problem l|7j is not trivial since it involves non-separable quadratic 
as well as non-smooth terms. Replacing {7} by the alternating mini- 
mization with respect to each vector leads to a variant of the so-called 
alternating direction method of multipliers (ADMM) 1111 : 

Algorithm ADMM 

1. Set k — 0, choose fii,fj,2 > 0, vo, wo, bo, and cq. 

2. repeat 

3. u fe+ i <- argmin^L||Gu-v fc -b fc ||| + ||u-w fe - c k \\? 2 



Vfc+l < 



M2 1 

argmin/ 2 (v) + f\ 

V 

- arg min/i (w) + 



Gu k+1 - v - bfc|| 2 

||w - Ufc+l 



Cfcllg 



bfc+i <— bfc + Guft + i — Vfc+i 
Cft+l <— c fe + u fe+ i — W fc+ i 
k <- k + 1 



9. until stopping criterion is satisfied. 

The proof of convergence in 1111 applies to a different variant 
of ADMM, which results from a different splitting. However, it is 
possible to show that this version can still be written as a standard 
ADMM and satisfies the conditions of the convergence theorem 1121 . 

3. PROPOSED METHOD 

3.1. Reformulation of the Problem 

The feasible set in problem l|2j is the ellipsoid 

B(e,B,y) = {x£l" : ||Bx - y|| 2 < e}, (10) 

possible infinite in some directions. Problem |2| can be written as 
an unconstrained (discontinuous) problem, 



(5) 



min 0(x) + Ie(6,i,o) (Bx - y) , 



(11) 



where is 



denotes the indicator function of set 5* C 
0, ifseS 



*-s(s) = 



+00, if s 4 S. 



(12) 



Notice that E(e, I, 0) is simply an e-radius Euclidean ball centered 
at the origin of R m . 

Since problem l | I I [ i clearly has the form {5}, its VS-based con- 
strained optimization reformulation is 

min <f>(w) + L E r Cs im(v), s. t. u = w (13) 

v = Bu y. 



3.2. Application of ADMM 

Performing the adequate translations (which are clear from compar- 
ing (6) with fl3}), the ADMM becomes the proposed C-SALSA. 

Algorithm C-SALSA 

1. Set k = 0, choose (11,(12 > 0, vo, wo, bo, and co. 

2. repeat 

3. u' <- w fc + c fc 

4. u" <- y + v fc + b fe 

5. u fc+ i ■ 

6- v' <- Bu fc+ i - y - b fc 

(e,I,0)l 

8. w'*-u fc+1 +c fc 

9. w fe+ i 

10. b fc+1 <- b fc + Bu fe+1 - y - v fc+1 

11. Cfc + 1 <- Cfc + U fc+ 1 - Wfc + 1 

12. k <- fc + 1 

13. until stopping criterion is satisfied. 



are min — II B u — \i"\\f, + II u — u' 

6 u M2 II 112 II 



'II-' 

2 



Vfc +1 <- argmin tj5( e ,i, ) (v) + ^||v - v'| 



'11 2 
2 



arg min </>(w) + ^ ||w — w 



'112 
2 



A key feature of C-SALSA is that the cost of each iteration is 
0(n log n), as confirmed by the following observations. Lines 3, 
4, 8, 11, and 12 simply involve adding vectors or scalars, thus have 
0(n) or O(l) cost. Line 5 consists in minimizing a strictly convex 
quadratic function, leading (with a = fii/^2) to 



Ufc+l 



«B T B 



I 



a B u 



+ u 



(14) 



As will be shown in Subsection |3.3[ in several cases of interest, this 
matrix inversion has 0(n log n) cost. Lines 6 and 10 involve matrix- 
vector products which, by the same reason, have 0{n log n) cost. 
Line 7 corresponds to the orthogonal projection of v' onto the e- 
radius £2 ball E(e, I, 0), which is an 0(n) operation: 



Vfc+i = V, 



Ev'/KHa, if 
if 



|x|| 2 > e, 
v|| 2 < e. 



(15) 



Finally, line 9 is simply Wfc+i = *0/n 2 ( w ') ( see QO- ^ <?H X ) = 
||x||i, the cost of Vl" is 0(n). If <j> is the TV norm, we use Cham- 
bolle's algorithm, which (although iterative) also has 0(n) cost [6|. 

3.3. Implementing (14) 

We will now show how (14) can be implemented with 0(n log n) 
cost in several cases of interest. If B represents a convolution, it is 
factorized as B = U H DU, where U is the unitary matrix (U s = 
U _1 ) representing the discrete Fourier transform (DFT) and D is a 
diagonal matrix. Thus, 



(aB T B + I)~ 



V" qIDI^ + I 



u, 



(16) 



where |D| 2 is the matrix with squared absolute values of the entries 
of D . Since a | D | 2 +1 is diagonal, its inversion costs O (n) . Products 
by U and U H have 0(n log n) cost, using the FFT algorithm. 



In frame-based regularization, the unknown image is represented 
on a frame (e.g., of wavelets or curvelets) and then the coefficients 
of this representation are estimated from the observed data, under 
some regularizer. A constrained formulation of this approach still 
has the form {2| but with different meanings for x and B: vector 
x now contains the frame coefficients of the unknown image Wx 
(the columns of W contain the elements of the adopted frame) and 
B = AW is now the product of an observation matrix A by the 
frame synthesis matrix W | 28 |. The only impact of this change on 
C-SALSA is in computing (14) , since AW is not diagonalizable by 
the DFT. This difficulty may be sidestepped under the assumption 
that W contains a 1-tight (Parseval) frame (i.e., W W H = I) and 
that A = U H DU, with D diagonal (e.g., a convolution). Using the 
matrix inversion lemma: 



aW ff A ff AW+I = I-W^A^IAA^ +I/a) AW 



(17) 

Since A = U fl DU, we have F = U H D* (|D| 2 + I/a) 1 DU, 
the computation of which has O(nlogn) cost, using the FFT to 
compute the products by U and U H . The cost of (17) will thus 
be either 0(n log n) or the cost of the products by W and W. For 
most tight frames used in image processing, there are fast 0(n log n) 
algorithms to compute these products 1211 . 

Finally, we considered the case of partial Fourier observations, 
which is used to model MRI acquisition and has been the focus of 
recent interest due to its connection to compressed sensing |5 |, 1 19 1. 
In this case, B = MU, where M is an m x n binary matrix (m < n) 
formed by a subset of rows of the identity, and U was defined above. 
Due to its particular structure, matrix M satisfies MM T = I; this 
fact together with the matrix inversion lemma leads to 



aB T B + 1 



-a/(l + a) U ff M T MU, 



(18) 



where M T M is equal to an identity with some zeros in the diagonal. 
Consequently, the cost of 1 18 1 is also 0(n logn). 



4. EXPERIMENTS 

All experiments were performed using MATLAB on a Windows XP 
laptop with a 2 GHz processor and 512 MB of RAM. 

We consider five standard image deconvolution benchmark prob- 
lems 1141 . summarized in Table [T| all on the well-known Camera- 
man image. We solve problem with (f>(x) = ||x||i (thus ^ is a 
soft threshold) and B = AW, where W is a redundant 4-level Haar 
wavelet frame and A is the blur operator. We set /ii — 112 and hand- 
tuned its value for fastest convergence. We compare C-SALSA with 
SPGL1 as follows. First, we run SPGL1 and then C-SALSA (from 
the same initialization), stopping when the constraint in (2) is satis- 
fied and the MSE of the estimate is below that obtained by SPGL1. 
Table [2] reports the number of iterations and CPU times taken in 
each of the experiments. Figure [4] plots the evolution of quadratic 
constraint || AWxfc — y|| 2 , in experiment 1. 

In the MRI reconstruction experiment, M models 22 radial ob- 
servations of the DFT and cj> is the TV norm [5|. Since SPGL1 can 
be used only for 0(x) = ||x||i, we compare C-SALSA with the code 
available in l\ -magic. Table [5] compares the 2 algorithms, in terms 
of computation time and the final MSE obtained. 



Table 1. Details of the image deblurring experiments. 



Experiment 


blur kernel 


<7 2 


1 


9x9 uniform 


0.56 2 


2A 


Gaussian 


2 


2B 


Gaussian 


8 


3A 


h IJ = l/(l + i 2 +j 2 ) 


2 


3B 


hij = 1/(1 + i 2 +i 2 ) 


8 



Table 2. Image deblurring using wavelets - Computation speed 



Experiment 


Iterations 


CPU time (seconds) 




SPGL1 


C-SALSA 


SPGL1 


C-SALSA 


1 


400 


136 


553.188 


118.953 


2A 


200 


152 


258.406 


130.203 


2B 


150 


120 


190.688 


115.375 


3A 


250 


57 


303.688 


48.5 


3B 


150 


46 


188.516 


40.5156 



5. CONCLUSIONS 

We have proposed a fast algorithm for solving constrained convex 
optimization problems usually known as basis pursuit denoising. 
Our algorithm is based on variable splitting and exploits augmented 
Lagrangian tools. Preliminary experiments with l\ and TV regular- 
ization show that the new algorithm outperforms existing methods in 
terms of computation time, by a considerable factor. Ongoing work 
includes a more thorough experimental evaluation of C-SALSA. 
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